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ABSTRACT 

We  describe  a  numerical  method  for  solving  the  steady,  three-dimensional,  incompress¬ 
ible  N'avier-Stokes  equations  in  cylindrical  geometry.  Also,  we  present  results  of  computa¬ 
tions  in  which  this  method  is  used  to  determine  the  flow  in  fluid-filled  cylinders  undergoing 
spinning  and  coning  motion.  Second-order  accurate  central  finite  difference  formulas  are 
used  to  approximate  derivatives  in  the  radial  and  axial  directions  and  a  Fourier  method 
is  used  to  approximate  the  angular  derivatives.  Nonuniform  grids  are  used  to  improve  the 
resolution  of  the  velocity  and  pressure  near  the  cylinder  walls.  The  system  of  difference 
equations  are  solved  using  an  iterative  method  based  on  successive-over-relaxation.  The 
method  has  been  found  to  be  very  efficient  in  terms  of  both  computer  time  and  storage. 
Results  of  the  numerical  method  applied  to  the  flow  in  spinning  and  coning  cylinders  are 
presented  for  several  cases  for  which  experimental  data  are  available.  In  addition,  pertur¬ 
bation  methods  are  used  to  study  the  data  at  small  coning  speeds  and  small  coning  angles. 
Numerical  results  of  this  no-coning  limit  are  compared  with  both  the  numerical  data  and 
experimental  data  at  low  coning  conditions. 
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SIGNIFICANCE  AND  EXPLANATION 


Numerical  methods  for  solving  for  the  flow  of  incompressible  fluids,  such  as  water, 
oil.  and  many  other  common  liquids,  are  important  for  the  study  of  many  engineering 
problems.  In  this  paper  we  describe  a  finite  difference  method  for  solving  the  equations 
of  steady  incompressible  flow  in  three  dimensions.  The  basic  method  can  be  applied  in 
any  coordinate  system,  but  since  the  application  to  be  studied  is  for  flow  in  a  cylinder, 
we  present  a  modified  method  which  takes  advantage  of  the  cylindrical  coordinate  system. 
The  numerical  method  is  used  to  solve  for  the  flow  in  fluid-filled  cylinders  undergoing 
both  a  spinning  and  coning  motion.  The  numerical  method  is  very  efficient,  running  on  a 
Vax  11  7*0.  even  though  the  flow  is  three-dimensional.  The  study  of  fluid-filled  cylinders 
undergoing  both  a  spinning  and  coning  motion  is  a  significant  problem  in  modern  ballistics, 
and  of  interest  to  L’.S.  Army  researchers  at  the  Ballistic  Research  Laboratory  and  Chemical 
Research  and  Development  Center.  Comparisons  of  the  results  of  computations  with  the 
new  scheme  with  experimental  data  show  the  scheme  to  be  very  accurate. 
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The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive  summary  lies 
with  MRC.  and  not  with  the  authors  of  this  report. 


A  NUMERICAL  METHOD  FOR  THE  INCOMPRESSIBLE 
NAVIER-STOKEJ  EQUATIONS  IN  THREE-DIMENSIONAL 
CYLINDRICAL  GEOMETRY 


John  C.  Strikwerda  A  and  Yvonne  M.  Nagel 

Introduction 

We  present  a  numerical  method  for  solving  the  steady,  incompressible  Navier-Stokes 
equations  in  three-dimensional  cylindrical  geometry.  The  difference  scheme  is  a  regularized 
central  difference  scheme  in  the  radial  and  axial  direction  as  introduced  by  Strikwerda  [17: 
with  a  Fourier  method  in  the  azimuthal  direction.  The  method  is  presented  as  applied  to 
the  computation  of  the  flow’  in  a  fluid-filled  cylinder  undergoing  both  a  spinning  and  coning 
motion.  The  values  of  the  velocity  components  and  pressure  are  assigned  to  a  common  grid, 
i.e.  a  non-staggered  grid  is  used,  and  grid  stretching  is  used  to  improve  the  resolution  near 
the  cylinder  walls.  Because  of  the  regularized  difference  method,  the  scheme  maintains 
second-order  accuracy  even  for  nonuniform  grids.  The  difference  equations  are  solved  by 
an  iterative  method  based  on  successive-over-relaxation  as  discussed  by  Strikwerda  jl8j 
and  thus  requires  only  one  three-dimensional  array  per  dependent  variable.  This  offers  a 
significant  savings  in  computer  storage  over  time-marching  methods.  The  iterative  method 
uses  line  successive-over-relaxation  together  with  an  exact  inversion  of  the  Fourier  operator 
for  each  line.  Each  line  consists  of  the  points  at  a  fixed  value  of  the  radial  and  axial 
coordinates. 

The  study  of  the  flow  in  fluid-filled  spinning  cylinders  is  of  importance  in  several  areas, 
especially  ballistics.  Stewartson  j  161  studied  the  inviscid  flow  in  a  spinning  cylinder  and 
Wedemeyer  22.23  applied  this  theory,  along  with  some  approximations,  to  the  case  of  high 
Reynolds  number  flows.  The  low  Reynolds  number  case  is  also  of  interest  since  projectiles 
filled  with  highly  viscous  fluids  have  exhibited  rapid  despin  and  growth  in  coning  angle.  In 
an  effort  to  understand  this  phenomenon  Miller  [9,10  and  D’Amico  and  coworkers  [2,14 j 
devised  experiments  in  which  fluid-filled  cylinders  were  subjected  to  simultaneous  spinning 
and  coning  motions.  Various  researchers  have  made  analytical  studies  of  the  problem,  e.g. 
Herbert  7.8  and  Murphy  11.12  .  A  numerical  study  was  done  using  a  finite  difference, 
time-marching  scheme  by  Vaughn.  Oberkampf.  and  Wolfe  21  using  a  method  of  Chorin 
3  .  The  present  method  was  applied  to  this  problem  because  of  the  need  for  a  more 
efficient  and  accurate  computer  code  for  studying  this  fluid  dynamics  problem. 

The  paper  is  organized  as  follows.  Section  1  describes  the  derivation  of  the  equations 
to  be  solved  by  the  numerical  method.  Section  2  presents  the  finite  difference  equations  and 
the  Fourier  method  for  the  incompressible  Navier-Stokes  equations.  The  iterative  solution 
algorithm  for  the  difference  equations  is  discussed  in  section  3.  A  perturbation  expansion 
in  terms  of  two  small  parameters  is  discussed  in  section  1  and  in  section  5  computational 
results  are  presented. 
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1.  The  Equations 

Consider  a  fluid-filled  cylinder  which  is  spinning  simultaneously  about  two  axes.  The 
first  axis  is  that  of  the  cylinder  and  the  second  axis,  the  coning  axis,  is  inclined  to  the 
cylinder  axis  by  a  fixed  angle  9.  Let  u.’  be  the  angular  velocity  with  which  the  cylinder 
rotates  about  the  coning  axis  and  let  1)  be  the  angular  velocity  of  rotation  about  the 
cylinder  axis.  The  fluid  velocity,  v  and  pressure.  P.  in  the  cylinder  are  governed  by  the 
incompressible  Navier-Stokes  equations  given  by 

—  =  -VP  ~  P_1  V2t~ 

Dt  (1.1) 

V  t7=  0, 

in  an  inertial  frame  in  which  the  coning  axis  is  fixed.  The  equations  and  variables  are 
non-dimensionalized  by  using  the  cylinder  radius,  a,  as  the  length  scale  and  the  product 
of  the  cylinder  radius  and  the  spin  rate  of  the  cylinder  about  its  axis  as  the  reference 
velocity,  hence  the  Reynolds  number  is  f la2  i/.  where  u  is  the  kinematic  viscosity  of  the 
fluid.  By  changing  to  the  non-inertial  coordinate  frame  which  rotates  with  angular  velocity 
jj  about  the  coning  axis  the  fluid  motion  becomes  steady.  It  is  also  convenient  to  compute 
the  difference  between  the  velocity  and  the  solid  body  motion  of  the  cylinder  rather  than 
the  velocity  itself.  Similarly,  it  is  convenient  to  compute  with  the  difference  between  the 
pressure  and  a  given  function  that  is  chosen  to  simplify  the  forcing  terms  arising  from  the 
coordinate  transformations. 

Based  on  these  transformations,  and  using  a  cylindrical  coordinate  system  aligned 
with  the  cylinder,  the  equations  describing  the  flow  are 


1  V2u  -r  (u-  V)u- 


J(t,  9)  x  u  +  Vp  =  2t  sin  9  r  cos  <i>  k 


V£=0  (1.3) 

where  r  is  the  ratio  of  the  angular  velocity  about  the  coning  axis  to  that  about  the  cylinder 
axis.  i.e.  jj  /  fl,  and  k  is  the  unit  vector  in  the  direction  of  the  z  axis,  the  cylinder  axis,  and 

J(t.  9)  —  2(-r  sin  9  cos  o.  r  sin  0sin  6. 1  +  t  cos  9)' 

=  2(1  ~  t  cos  9) k  -  2 t  sin  9  i{<t>) 

in  the  cylindrical  coordinate  system,  where  i (<p)  is  the  unit  vector  in  the  i  direction.  The 
j  s  plane  is  the  plane  of  the  two  axes  which  is  also  the  plane  <p  --  0.  The  velocity  u  in 
equations  (1.2)  and  (1.3)  is  relative  to  the  solid  body  rotation  of  the  cylinder  and  in  the 
coordinate  system  rotating  with  the  coning  motion,  thus  the  boundary  condition  for  the 
system  (1.2.  1 .3)  is 

v  =  0  (1.4) 

on  the  cylinder  boundary.  The  relations  between  the  actual  velocity  and  pressure  (v.P)  in 
equation  (1.1)  and  the  computed  velocity  and  pressure  (u.p)  in  equations  (1.2)  and  (1.3) 


are 


r 


(1.5) 


V  P 


(1.6) 


—  - r 2  -  r2r  cos  6  —  —  ((r  cos  0 cos  6  z  sin  0)2  +  r2sin2#), 


where  t  is  the  vector  in  the  direction  of  the  coning  axis  with  magnitude  r.  The  system 
(1.2,  1.3)  holds  for  0  <  r  <  1  and  -b  <  z  <  b  where  6  is  the  aspect  ratio  of  the  cylinder, 
defined  by  the  ratio  of  the  length  of  the  cylinder  to  its  diameter. 

In  addition  to  the  system  (1.2.  1.3).  it  is  also  of  interest  to  examine  the  equations 
resulting  from  a  linearization  about  r  and  B  equal  to  zero,  as  will  be  discussed  in  section 
5.  The  resulting  system  is 


-R~lV2u 


dH 

- 2fcx«  +  Vp=2r  cos  0  k 

d0 


(1.7) 


Vu  =  0.  (1.8) 

This  system  will  also  be  used  to  explain  the  iterative  solution  algorithm  which  uses  a 
linearization  to  determine  the  updates  to  the  solution. 


2.  The  Numerical  Method 

The  numerical  method  to  solve  the  equations  of  flow  is  based  on  the  regularized  finite 
difference  method  of  Strikwerda  jl7j  together  with  a  Fourier  or  pseudo-spectral  method. 
Finite  differences  are  used  to  approximate  derivatives  in  the  axial  and  radial  direction  and 
the  Fourier  method  is  used  to  approximate  derivatives  in  the  azimuthal  direction.  The 
Fourier  method  is  much  more  accurate  than  the  finite  difference  method  for  periodic  in¬ 
dependent  variables  such  as  0,  (Gottlieb  and  Orszag  5;),  and  this  allows  for  a  substantial 
savings  in  computer  storage  and  thus  also  in  computer  run  time.  The  regularized  differ¬ 
ences  make  it  possible  to  compute  both  the  velocity  and  pressure  on  the  same  grid,  as 
opposed  to  staggered  grid  methods.  Also,  grid  stretching  is  used  to  place  more  points 
in  the  regions  near  the  container  walls  to  better  resolve  the  flow  field.  The  regularized 
differences  allow  this  to  be  done  simply  without  loss  of  accuracy.  The  nonuniform  grid  is 
defined  with  the  use  of  mappings 


r  =  ap-  ( 1  -  a)p3  (2.1) 

*  =  6(0f  +  (l  -5)(l5f3-8f5)/7),  (2.2) 

which  map  the  region  0  <  p  <  1  and  —  1  <  c  <  1  onto  the  cylinder.  Because  the  resolution 
at  the  end  walls  was  critical  to  the  accuracy  of  the  solution  the  stretching  (2.2)  was  chosen 
to  have  the  property  that  d2  z> d<;2  vanishes  at  <  equal  to  il  so  as  to  obtain  a  more  uniform 
stretching  there.  The  values  of  a  and  3  are  chosen  to  make  ar!dp  and  dz/d<;  less  than  1 
and  6.  respectively,  at  p  and  C  equal  to  1.  The  variables  p,  C-  and  the  angular  variable,  o 
are  discretized  using  uniform  grid  spacing,  i.e. 

A P=  !/(/-!)•  P,  =  (i  -  1) Ap, 

Ac  --  2 f(K  -  1).  sk  =  {k-  1 ) Ac. 

A<Z>  =  2n  J,  Oj  -  ( j  -  1)A®. 


o 

«> 


The  differentia)  equations  were  written  in  terms  of  the  new  independent  variables  p  and  f 
and  then  differenced  using  central  difference  formulas.  As  an  example, 


1  d_ 
r  dr 


I  du  \ 


1  d  r  du 
rr'(p)  dp  r’(p)  dp 


(2.3) 


1  /r«^l  2,  ,  r»  — 1,2  /  ,  \  1 

{ui^\,).k  Ui.j.k)  ,  (ut,j.k  U,  -  l.;./c))  .  ,  • 

r,ri  ,-M/2  1  —  1/2  “^P 


wher^  r,  =  r(p,).  r'  =  r'(pt),  etc.  The  gradient  terms  and  the  divergence  terms  were 
differenced  using  regularized  differences  as  introduced  by  Strikwerda  17  .  For  example, 
dp  dz  in  (1.2)  and  du  'dr  in  (1.3)  were  approximated  by 


<^P  ^  2  jPijA4-!  Pi.j.k-l  Pi.j.k+2  ~  3pttJ|fc+ 1  Pt,j,/c-l, 

~  s'(f*j  2Af  6A<; 


dru 

rdr  ,,j> 


1  /  I  U«+  1  ,j,k  ri  -  1  tt  i  —  l  tjtk  \  3u,  J.lt  ■+-  3u,_  1,;,*  Ut  —  2,],k 


rr'(Pi 


r- 


2Ap 


)- 


6r'(p,)Ap 


As  discussed  in  Strikwerda  17  ,  these  regularized  differences  were  used  only  on  the  terms 
dpi  dr  and  dpidz  in  the  pressure  gradient  and  on  the  terms  dru!  rdr  and  dwidz  in  the 
divergence  equation.  The  approximations  to  dp/dcp  and  dvjd<t>  will  be  discussed  later. 

In  the  Fourier  method  approximations  to  derivatives  with  respect  to  the  angular  vari¬ 
able.  0.  are  obtained  as  follows,  using  the  pressure  as  an  example.  For  each  choice  of  the 
radial  and  axial  grid  indices,  i  and  k.  the  discrete  function  p,,;>  is  represented  using  the 
discrete  Fourier  transform  as 


J  2 

P.,;.*  =  51  °!Jt).fcsinn*J  -  bun.kcosnQj-  ( 2  4 ) 

n-0 

where  aj*0  k  and  ot4J/2  *.  are  an^  the  prime  on  the  summation  indicates  that  the  first 

and  last  terms  are  weighted  with  a  Similar  expressions  with  coefficients  k 

hold  for  the  velocity  components  u.v.  and  w.  for  l  -  1.2.3.  respectively.  Considering  the 
right-hand  side  of  (2.4)  as  a  continuous  function  of  0.  and  then  differentiating  this  function 
with  respect  to  <z>  we  obtain 

J,  2 

do  *  Y,  no!4n.*cosn<Zb  -  nbl4lksmno}.  (2.5) 

'•J-k  n->J 

I  he  coefficients  and  b{4^  k  are  easily  obtained  by  the  use  of  formulas 


a 


(41 

i  ,n,k 


2 

J 


j  i 

5^  p*,j*  sin  n<j>j 


4 


and 


» (4)  2  ^ 

6.,n>  =  j  2s  pt,hk  COS  n<i>j. 

]=<' 

To  maintain  the  regularity  of  the  difference  method,  (see  Bube  and  Strikwerda  1  ) 
the  approximations  to  dp  d<t>  in  the  pressure  gradient  in  ( 1 .2)  and  dv ; d<p  in  the  divergence 
equation  (1.3)  were  approximated  as  in  (2.5)  with  the  addition  of  the  term 

(2-6) 

to  (2.5)  and  the  subtraction  of  a  similar  term  to  compute  dv/d<t>.  The  use  of  the  factor 
!  was  somewhat  arbitrary,  however  the  omission  of  this  term  (2.6)  led  to  nonconvergence 
of  the  solution  in  almost  all  cases.  The  importance  of  such  a  term  can  be  seen  in  that 
without  the  term  (2.6)  the  Fourier  mode  b, tj/2.k  cos  o}  which  is 

will  not  contribute  to  the  derivative  (2.5).  Thus  the  scheme  would  allow  grid  scale  oscilla¬ 
tions  in  the  <t>  direction.  The  use  of  this  extra  term  in  occurences  of  first  derivatives  of  the 
velocity  components  with  respect  to  <j>  in  (1.2)  or  (1.7)  resulted  in  poorer  performance  of 
the  iterative  solution  algorithm  and  so  was  not  used.  The  effect  of  the  regularizing  terms 
has  only  to  do  with  removing  high  frequency  oscillations  and  does  not  affect  the  formal 
accuracy  of  the  method. 

Along  the  z  axis  of  the  cylindrical  coordinate  system  the  values  of  the  velocity  and 
pressure  were  determined  by  interpolation  from  the  neighboring  grid  points.  In  particular, 

u'lj.tc  =  (4«?2,fc  -  U>3,k)'3 


and 

where 


and  similarly  for  p,  *.  Note  that  Wij.k  and  P],},k  are  independent  of  j.  Because  of  the 
multiple-valued  representation  of  the  velocity  along  the  axis  of  the  cylindrical  coordinate 
system  we  have 

«t,;,tc  =  t/i>COS0;  -  Vj.fcSinOj 
tlj./c  — Il.fr  COS  O  j 


where  U\k  and  V\ are  the  velocity  components  on  the  axis  in  the  cartesian  coordinate 
system.  The  values  of  b'l,*:  and  Vi^  are  given  by 


U  \  ,k  ~  (4t'2,tc  ~  f  3.(c)  3 

l'l .*  “  (4l'2.*  “  ^3.k)  3 


o 


with 


f,  1  XT' 

c  i.k  —  /  Ut ,j.fc  cos  0j  —  vt  j,t  sin 


and 


1 

l  t.fc  =  j  >  ^..j.itsin^  -  w.-j.fc  cos  </>_,. 


j  =  i 


As  shown  in  Strikwerda  and  Nagel  20  this  treatment  at  the  origin  preserves  the  second- 
order  accuracy  of  the  method. 

The  values  of  the  pressure  on  the  cylinder  walls  were  determined  by  extrapolation 
from  the  interior  by  the  formulas 


Pl,j,k  =  3pi-l,]  k  -  ^Pl-2,}.k  ~  Pi -3,j,k 


and 


Pi,j,K  -  3p,j,K_i  -  3plyJ.K-2  -  Pi,j ,K  -  3 


and  similarly  for  p1}  j. 

3.  The  Iterative  Solution  Algorithm 


The  system  of  nonlinear  difference  equations  is  solved  using  a  modified  line  successive- 
over-relaxation  method  (LSOR).  The  basic  method  is  described  in  Strikwerda  18!  for  the 
case  with  linear  finite  difference  equations.  Because  the  Fourier  method  is  used  in  the 
azimuthal  direction  it  is  most  convenient  to  use  line  SOR.  with  each  line  being  the  grid 
points  in  the  azimuthal  direction  for  each  given  value  of  the  radial  and  axial  coordinates. 
The  coupling  between  the  three  velocity  components  and  the  natural  periodicity  of  the 
azimuthal  coordinate  leads  to  a  periodic  system  of  equations. 

To  describe  the  solution  algorithm  it  is  best  to  regard  the  discrete  solution 


(Ut,J.k-  V,  U\  j_k.  P,.j.k) 

as  continuous  functions  of  o  with  a  finite  Fourier  expansion,  that  is.  as 


(Ut,fc(0),  »',,*(©).  U.'t,Jc(0).p, >(<?))  --  (u, ,*(<?),  p,. *(<£)). 


The  finite  difference  and  Fourier  scheme  for  the  linearized  problem  (1.7)  can  be  written  as 


1  .k 


xi.K  -  4 t,k  -  4t  .k  -  4t,/c  - 

ut  l.k  ~  \.kU^  1  .*  "  Alk  ]ut.k  I  ~  -v  I  u  +  1  A^kUi.k 


+  2.kPi~2.k  ~  Gt* \,kPi~  l.*  *  67  ,  kP,  l.k  d  ^  + 2^' .*  + 


ii,k 


(3.1) 


~  6');7,P,.*  1  "  G\-*kpuk  =  2r,  cos  <f>}k 
(> 


-»t  .k 


i  .k 


Jl 


where  each  of  A\'kk  and  G\'k  are  matrices  of  differential  operators  in  <J>.  In  particular,  the 
operator  A['k  is 


where  L,  *  is  the  operator 


(  R~l{L<k  -R-'bl-*  o 

-R  -2  o 
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The  iterative  method  is  then  written  as 


A  l  ,K  /  —  J/  -t-  1  -*1/  \  - ► 

res>,* 


(3.2) 


where  res,,*  is  the  difference  of  the  right-hand  and  the  left-hand  sides  of  equation  (3.1) 
evaluated  using  the  most  current  values  of  u,  and  u  is  the  SOR  iteration  parameter.  The 
equation  (3.2)  was  solved  by  computing  the  Fourier  coefficients  of  the  components  of  the 
velocity  residual  vector  res,  *  for  the  current  value  of  (i.k)  and  solving  for  the  Fourier 
coefficients  of  the  velocity  update  ul'~k  1  -  ti”  k.  This  leads  to  a  set  of  linear  equations  for 
the  Fourier  coefficients  of  the  velocity  update  which  are  solved  by  Gaussian  elimination. 
Each  linear  system  is  of  size  at  most  six  by  six  for  the  sine  and  cosine  terms  for  each  value 
of  n  as  in  (2.4)  for  the  three  components  of  the  velocity  update. 

For  the  nonlinear  system  of  equations  the  same  algorithm  was  used  with  res,  *  being 
computed  with  the  addition  of  the  nonlinear  convection  terms. 

After  the  velocity  was  updated  at  all  the  grid  points  by  one  pass  of  the  LSOR  operator, 
the  pressure  was  updated  bv  setting 


h  ■  u 


- 1 


(3.3) 


w'here  is  an  iteration  parameter  as  described  in  Strikwerda  18  and  V*-  is  the  discrete 
divergence  operator.  As  in  Strikwerda  17,18  .  we  did  not  attempt  to  satisfy  the  equation 


V h  ■  =  0.  (3.4) 

but  rather 

^ h  ■  Ut,],k  ~  (3.5) 

where  6*  is  the  average  value  of  the  discrete  velocity  divergence.  In  all  the  cases  presented 
here  the  value  of  bh  w.is  on  the  order  of  the  truncation  error.  As  discussed  in  Strikwerda 
17  ,  the  difference  equation  (3.4)  with  the  boundary  conditions  (1.4)  will  in  general  be 
an  over-specified  system  and  need  not  have  a  solution.  By  allowing  the  divergence  of  the 
velocity  field  to  be  a  constant,  but  nonzero  value,  the  system  has  a  unique  solution. 


This  method  has  the  advantage  that  only  one  three-dimensional  array  is  needed  for 
each  of  the  velocity  and  pressure  variables  as  opposed  to  time-marching  methods  which 
require  two  such  arrays.  The  LSOR  iteration  has  a  good  rate  of  convergence,  allowing  for 
the  solutions  to  be  computed  on  a  VAX  11  '780  computer  in  at  most  several  hours  of  cpu 
time. 

4.  Perturbation  Expansions  in  the  Coning  Speed  and  Angle. 

Many  of  the  experimental  results  are  for  small  values  of  the  parameters  r  and  0, 
therefore  it  is  useful  to  make  a  perturbation  expansion  of  the  solution  of  the  system  (1.2) 
and  (1.3)  in  terms  of  these  parameters.  Since  the  forcing  term  is  proportional  to  t  sin#  the 
velocity  and  pressure  are  also  proportional  to  this  quantity.  Consider  then  the  expansion 
in  the  form 

(u,  p)  —  t  sin 6  Y,  (urn  n,pm'n)Tms\r\n  6.  (4.1) 

m.n  =  0 

For  each  term  of  the  expansion  we  have 


V  •  un'm  =  0. 

and  ach  un-m  vanishes  on  the  boundary.  The  equation  for  (u°’°.  p0,0)  =  (t7°,p°)  is  the 
same  as  (1.7) 

di/-° 

-R  ‘V2u°  *  — - 2k  x  u°  -r  Vpu  -  2rcos<j>  k.  (4.2) 

OQ 

This  equation  has  been  analyzed  by  Gerber  et  al.  4  for  the  case  of  high  Reynolds  number 
using  separation  of  variables  and  numerical  methods.  If  we  let  the  operator  on  left  side  of 
equation  (4.2)  be  denoted  by  L  and  set  Wm  n  =  (um’n.  pm'n).  we  obtain 

L\Yi  0  ~  -2k  *  v:'  (4.3) 

and 

IM’ul  =  0 

hence  It'01  -  0. 

Since  the  operator  L  is  linear  and  the  forcing  term  in  (4.2)  involves  only  the  k  =  1 
Fourier  mode  we  see  that  IV  and  W1"  contain  only  the  k  —  1  Fourier  mode.  Similarly 
lor  IV  1  2  and  IV2"  which  satisfy 

LW  1  2  -  k  >:  u  \  LW 2  ,1  =  -2k  *  u1'0. 

The  equation  for  IV  1  1  has  a  right-hand-side  involving  a  quadratic  in  u0  and  the  vector 
i(o)  times  if  .  and  thus  W  1,1  contains  only  the  Fourier  modes  for  k  =  0  and  2.  Similarly, 
11  21  also  has  only  these  modes. 

Nusca  et  al.  14)  defined  the  coefficient  of  pressure.  C,,.  as 

max  P]  (r.  z)  6  ~  1 
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for  values  of  ( r.z )  on  the  end  wall  of  spinning  and  coning  cylinders,  where  pi(r.z)  is  the 
amplitude  of  the  Fourier  mode  for  k  -  1.  That  is,  P\{r.z )  is  determined  by 


p(r,0.z)  =  ^  pk(r,z)s\n(ko  -  (*)■ 


k=o 


from  the  perturbation  expansion  in  r  and  sin#,  and  considering  also  the  solid  body  con¬ 
tribution  to  the  true  pressure,  e.g.  (1.6),  we  obtain 


C,,(r.z)  ^  T{(p\f(r,z)  -  7p]'^’(r,  z))2  -  (Pi,c(r,z)  +  TP\\Jc[r,z)  -  ttz )2) 


k  2  v  1/2 


0(r2  sin2  #,  r3). 


(4.4) 


Note  that  the  terms  for  r2  sin  #  arising  from  W  1,1  do  not  affect  Cr  since  W  1,1  involves  only 
the  Fourier  modes  for  k  -  0  and  2.  Also  we  set  Cv  to  have  the  same  sign  as  r. 

The  coefficients  in  (4.4)  are  easily  computable.  The  coefficients  p\  s  and  p,  are 
computed  by  solving  the  equations  (4.2).  The  coefficients  p\'^  and  pj can  be  obtained 
either  by  solving  the  system  (4.3)  or.  as  was  done  here,  by  solving  the  system 


q  -0 

-  /?”  1  V  2  -  —  0/1  - l/T  .  .v0  r?_f' 

do 


2(1  -  r)fc  x  u  '  -  Vpl  -  2rcos<p  k. 


(4.5) 


This  system  is  obtained  by  performing  an  expansion  of  the  solution  in  terms  of  sin#  only, 
similar  to  the  way  (4.2)  was  obtained.  The  system  (4.5)  was  solved  for  #  equal  to  0.0 
and  for  r  equal  to  both  a  small  positive  and  a  small  negative  value.  Accurate  estimates 
of  p }  ‘ and  p } can  then  be  obtained  by  computing  the  divided  difference  of  the  Fourier 
coefficients  of  the  pressure  on  the  endwall  using  both  values  of  r. 

The  coefficients  of  liquid  roll  moment,  C/rm.  and  the  coefficient  of  liquid  side  moment, 
•  see  Murphy  11.12  ,  can  also  be  developed  as  expansions  in  r  and  sin#  using  the 
expansion  (4.1).  The  expansions  for  C;rrn  and  are 

C,rm  =  rTj  •'(/?,  6)/ 2tr6  +  0(r2) 


and 

C,.m  =  ri0’°(/?.#)/27r6-O(r), 

where  Tj  1  and  Tj'"  are  terms  in  the  expansions  of  the  moments  of  force  resulting  from 

(4.1) .  The  yaw  moment  T't ''"(/?. 6),  and  thus  C^m.  is  easily  obtained  by  solving  the  system 

(4.2) .  The  despin  moment  T'  ]  is  more  difficult  to  obtain  and  no  attempt  was  made  to 
compute  it. 

5.  Computational  Results 

The  numerical  method  described  above  has  been  used  to  compute  the  flow  in  a  spin¬ 
ning  and  coning  cylinder  for  a  large  number  of  cases.  Here  we  present  the  results  of  a 


9 


few  representative  runs.  A  discussion  of  other  cases  and  their  engineering  significance  will 
appear  elsewhere. 

Figure  1  displays  the  data  for  Cp  measured  on  the  endwall  of  the  cylinder  for  the  case 
R  -  3.1,6  =  3.148  and  6  =  2°.  Figure  la  shows  Cp(r,  6)  at  the  radial  location  r  =  0.667 
and  Figure  lb  shows  the  same  quantity  at  r  =  0.434.  The  results  of  the  computations  are 
marked  with  a  x  and  the  gyroscope  data  of  Hepner  !6l  are  marked  with  a  o.  In  addition 
to  the  data  points,  the  curve  based  on  the  expansion  (4.4)  for  Cp  is  also  displayed.  For 
positive  values  of  t  there  is  excellent  agreement  between  all  three  sets  of  points,  however, 
for  negative  values  of  r  the  experimental  data  differs  from  the  computed  data  and  the 
asymptotic  formula.  Note  that  the  computed  solution  and  the  curve  based  on  (4.4)  are 
not  completely  independent  since  the  numerical  method  to  compute  the  coefficients  in 
(4.4)  is  essentially  the  same  as  the  basic  numerical  method.  The  significant  qualitative 
difference  between  the  experimental  data  and  the  other  data  for  negative  values  of  r  may 
indicate  a  systematic  error  :  'l  .perimental  data.  Analytical  results  of  Sedney  (15] 
agree  qualitatively  with  computational  results. 

Figure  2  displays  the  coefficient  of  pressure  Cp  on  the  endwall  as  a  function  of  radius 
for  0  equal  to  2  degrees  and  two  different  values  of  r.  Figure  2a  shows  results  for  r  =  0.350 
and  Figure  2b  displays  results  for  r  =  0.100,  the  experimental  data  and  numerical  results 
are  marked  as  in  Figure  1.  In  both  figures  the  excellent  agreement  between  the  basic 
numerical  results  and  the  perturbation  results  verify  the  essential  linearity  of  the  problem 
for  small  values  of  r  and  6. 

These  calculations  were  made  with  7  -  11,  J  —  6,  K  —  33,  with  grid  stretching 
parameters  a  and  3  chosen  so  that  the  dr /dp  was  0.8  at  r  =  1.0  and  dz/di '  was  0.86  at 
z  -  b.  The  run  with  r  =  0.350  took  82  iterations  until  the  f2  norm  of  the  changes  in  each 
variable  was  less  than  10-4.  This  took  approximately  550  seconds  of  cpu  time  on  the  Vax 
11/780  at  the  Mathematics  Research  Center. 

Figure  3  displays  Cjgm  as  a  function  of  the  Reynolds  number,  R.  for  the  case  6  =  2.20 
for  R  between  0  and  350.  Since  Cjsm  is  related  to  the  yaw  moment,  a  negative  value  of 
this  quantity  indicates  a  tendency  to  decrease  the  angular  velocity  about  the  coning  axis 
and  a  positve  value  indicates  a  tendency  to  increase  this  angular  velocity.  The  variation 
in  the  sign  of  Cism  with  increasing  R  is  due  to  the  reduction  of  the  viscous  effects  on 
the  moment  as  R  is  increased.  Indeed,  the  contribution  to  of  the  pressure  on  the 

side  wall  is  of  positive  sign,  the  other  contributions,  i.e.  pressure  on  the  end  wall  and  the 
viscous  contribution,  are  all  negative  for  this  value  of  the  aspect  ratio.  As  the  Reynolds 
number  increases  the  viscous  side  wall  contribution  decreases  in  magnitude  changing  the 
sign  of  Cfsm .  A  significant  factor  in  the  flight  instability  of  liquid  filled  projectiles  is  an 
increase  in  the  coning  angular  velocity  and  is  most  likely  related  to  the  positive  values  of 
C®gm  shown  in  Figure  3. 

6.  Conclusions 

The  numerical  method  presented  in  this  paper  has  been  shown  to  be  an  efficient 
and  accurate  method  for  computing  solutions  to  the  steady  incompressible  Navier-Stokes 
equations  in  three  dimensions.  The  method  has  been  applied  to  the  computation  of  flows  in 
cylinders  undergoing  spinning  and  coning  motion  and  results  agree  well  with  experimental 
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Cp  vs.  r  at  radius  =0.667 


Tau 

Figure  la 

Cp  vs.  r  at  radius  =0.434 


Tau 

Figure  lb 


n 


data.  The  use  of  the  perturbation  analysis  corroborates  the  accuracy  of  the  calculations. 
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